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Abstract 


Electrodynamic  Tethers  (EDT’s)  offer  a  real  option  for  zero-propellant  orbital  maneuvers  in 
the  near  future.  By  controlling  the  electrical  current  through  a  long  conductive  cable  aligned  with 
the  local  vertical  and  in  the  presence  of  a  magnetic  field,  the  tether  experiences  an 
electrodynamic  thrust.  The  local  ionosphere  provides  the  necessary  electrons  for  the  generation 
of  an  electrical  current.  Previous  investigation  has  been  focused  on  feed-forward  or  open-loop 
control  schemes.  Open-loop  control  methods  are  very  susceptible  to  model  error.  The  relevant 
models  for  an  EDT  system  are  the  atmospheric  density  model  and  the  magnetic  field  model.  This 
paper  will  be  concerned  with  errors  in  the  atmospheric  density  model.  The  problem  consists  of 
two  parts:  solving  the  open  loop  non-linear  optimal  control  problem  and  solving  the  associated 
linear  feedback  system  to  generate  a  control  law. 

To  solve  the  first  part  we  assume  the  orbit  remains  nearly  circular.  We  apply  the  method  of 
averaging  to  the  state  dynamics  to  track  secular  changes  only.  The  short  period  motion  of  the 
spacecraft  drives  the  shape  of  the  control.  We  vary  the  coefficients  on  a  five  term  modified 
Fourier  series  describing  the  tether’s  alternating  electrical  control  current.  The  series  is  modified 
by  using  square  waves  rather  than  sine  and  cosine  waves.  The  open-loop  control  solution  is  then 
used  as  reference  in  the  feedback  problem. 

Solving  the  associated  linear  feedback  system  involves  linearizing  the  state  dynamics. 
Standard  linearization  yields  a  classic  state-space  structured  system  using  state  error  to  generate 
control  corrections.  We  assume  complete  state  feedback  to  simplify  the  solution.  Treating  the 
system  as  linear  time  invariant,  we  update  the  gain  matrix  once  per  orbit. 

Results  indicate  this  strategy  improves  performance  and  reliability  of  a  system  with  model 
errors  and  un-modeled  disturbances,  particularly  for  maneuvers  that  remain  in  the  LEO  regime 
for  an  extended  time. 


Keywords:  EDT,  Electrodynamic  Tether,  Optimal  Control,  Feedback  Control,  model  error 
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An  Approach  to  Optimal  Control  of  Electrodynamic  Tethers 
in  a  Stochastically  Varying  Drag  Environment 

Alexander  J.  Buck 

I.  Introduction 

The  Low  Earth  Orbit  (LEO)  environment  provides  many  scenarios  where  maneuvering  in 
orbit — perhaps  for  an  extended  duration — is  very  necessary.  The  tenuous  atmosphere  is 
constantly  slowing  down  spacecraft  through  drag,  lowering  their  altitude;  orbital  debris  can 
travel  at  relative  speeds  in  excess  of  10  km/s,  fast  enough  to  damage  or  destroy  a  spacecraft;  and 
spacecraft  are  constantly  being  perturbed  out  of  their  orbit  by  small  but  constant  effects  such  as 
solar  wind,  radiation  pressure,  slight  variations  in  the  gravitational  field,  and  the  gravitational 
influence  from  the  Sun  or  Moon.  For  example,  the  International  Space  Station  loses  about  90 
meters  of  altitude  per  day  due  to  drag  [1],  At  that  rate  Station  would  lose  over  30  km  in  one  year 
which  is  10%  of  its  orbit  altitude.  To  prevent  the  station  from  re-entering  the  atmosphere  and 
burning  up  after  just  a  few  years,  every  few  months  resupply  modules  are  launched  with  extra 
fuel  to  boost  the  station  back  up  to  a  higher  altitude.  Orbital  debris  removal  is  another  mission 
requiring  large  amounts  of  orbital  maneuvering  for  extended  periods  of  time.  In  early  2009  the 
Cosmos  225 1 — a  defunct  Russian  communication  satellite — collided  with  and  destroyed  Iridium 
33 — an  operational  US  communications  satellite.  This  event  added  roughly  1300  new  pieces  of 
debris  larger  than  roughly  5cm  radius — the  approximate  detectable  radius  of  the  US  Space 
Surveillance  Network  sensors  [2].  About  once  per  year  the  ISS  has  to  perform  small  maneuvers 
to  avoid  this  debris,  as  it  did  in  October,  2010  when  a  piece  from  a  decommissioned  NASA 
satellite  broke  off  and  had  a  close  approach  with  the  station.  A  0.4  m/s  boost  maneuver  was 
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executed  to  move  the  station  away  from  the  debris  [3].  Instead  of  avoiding  space  debris,  there  are 
proposed  missions  to  actively  rendezvous  with  space  debris  in  order  to  clean  up  the  LEO 
environment.  This  mission  profile  has  a  very  high  maneuverability  requirement  in  order  to 
rendezvous  with  and  capture  debris  in  all  its  various  orbits.  These  three  reasons  only  begin  to 
illustrate  why  we  need  to  have  maneuvering  capabilities  in  orbit. 

Chemical  propulsion  systems  are  currently  the  workhorse  for  this  capability.  These  systems 
are  largely  reliable  but  there  are  certain  problems  with  the  technology.  Chemical  rockets  rely  on 
a  reaction  mass  to  produce  the  momentum  exchange  and  pressure  difference  that  generate  thrust 
but  are  limited  by  the  amount  of  propellant  they  can  carry.  Due  to  the  disturbances  previously 
mentioned,  a  satellite  with  no  maneuverability  becomes  a  large  piece  of  orbital  debris,  or  a 
“zombiesat”  as  was  clearly  shown  in  April,  2010  when  Galaxy  15  turned  unresponsive  and 
began  drifting  towards  its  neighboring  satellites  [4],  Therefore  before  a  satellite  loses  the  ability 
to  maneuver  it  is  disposed  of  by  re-entry  or  placement  into  a  graveyard  orbit.  The  more  a 
spacecraft  has  to  maneuver  the  faster  it  uses  its  available  propellant,  the  shorter  its  lifespan  is. 

This  can  be  problematic  because  satellites  can  be  large  and  expensive  and  we  want  them  to  be 
operational  for  a  very  long  time,  often  15  years  or  more.  To  support  a  long  mission  life,  a 
significant  amount  of  propellant  (mass)  must  be  launched  with  the  satellite.  This  extra  mass  must 
be  launched  and  stored  in  the  spacecraft.  Considering  that  it  costs  roughly  $20,000  [5]  to  put  one 
kilogram  of  mass  into  orbit,  each  kilogram  of  propellant  can  add  significantly  to  launch  costs. 

Chemical  rockets  have  been  around  for  a  long  time  and  their  fuel  efficiencies — measured  by 
specific  impulse  (Isp) — have  gotten  better  over  time,  but  physical  limitations  of  the  combustion 
chemistry  have  more  or  less  limited  the  best  chemical  motors  to  a  vacuum  specific  impulse  of 
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455  seconds.  The  specific  impulse  of  a  motor  is  directly  linked  to  how  much  velocity  change  it 
can  induce  via  the  ideal  rocket  equation, 

&V  =  Ispg  In  (»!,/»!,) 

The  subscripts  i  and  /  denoting  initial  and  final  masses,  the  difference  being  the  propellant 
mass.  For  a  mission  that  requires  a  total  orbital  velocity  change  (AV)  of  10  km/s  with  a  vacuum 
Isp=455  seconds,  90%  of  the  initial  launch  mass  must  be  fuel.  This  means  a  100  kg  spacecraft 
requires  900  kg  of  fuel  to  be  launched  along  with  it  [6].  Using  the  $20,000/kg  mark  cited  earlier, 
the  total  launch  cost  for  this  system  is  roughly  $20  million. 

There  are  other  technologies  being  developed  and  used  to  help  address  some  of  the 
complications  and  drawbacks  of  chemical  propulsion  systems.  One  such  technology  is  electric 
propulsion,  which  accelerates  heavy  ions  through  an  electromagnetic  field  to  produce  thrust. 
Electric  propulsion  still  depends  on  a  reaction  mass  to  generate  thrust.  This  means  there  is  still  a 
finite  fuel  life,  albeit  considerably  longer  given  the  orders  of  magnitude  higher  Isp  values  for 
these  systems — on  the  order  of  5,000-10,000+  seconds.  For  the  same  10  km/s  AV  mission 
requirement  an  electric  propulsion  motor  operating  at  5,000  seconds  Isp  can  achieve  that  with 
only  20%  of  the  initial  payload  launch  mass  as  propellant.  To  operate  a  100  kg  spacecraft  would 
require  25kg  of  propellant  to  be  launched  along  with  it  reducing  the  total  launch  mass  by  a  factor 
of  8.  As  good  as  electric  propulsion  is  there  is  another  propulsion  technology  being  tested  called 


electrodynamic  tethers  (EDTs). 
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Figure  2  -  An  Electrodynamic 
Tether  (EDT).  Credit  Naval  Research 
Labs 


Electrodynamic  tethers  offer  the  real  prospect  of  zero- 
propellant  maneuvering  in  the  earth  orbit  environment.  This  is 
achieved  by  using  a  very  long  conductive  tether  (Figure  2) — hence 
the  name  electrodynamic  tether — to  interact  with  the  geomagnetic 
field  via  moving  electrons  (current).  The  resultant  force  on  the 
conductive  tether  as  described  by  the  Lorentz  force,  the  force  on  a 
current  carrying  wire  in  a  magnetic  field,  is  a  propulsive  force.  The 
Lorentz  force  is  defined  as  F  =  7L.rB  where  B  is  the  local  magnetic 


field,  L  is  a  vector  along  the  length  of  the  tether  and  I  is  the  control  current  causing  the  force.  In 
Figure  1  the  force  is  directed  out  of  the  page.  This  is  very  different  from  the  momentum 
exchange  through  particle  acceleration  that 
both  chemical  and  electric  propulsion  thruster 
designs  operate  under.  Similar  to  electric 
propulsion  thrusters,  EDTs  can  only  generate  a 
very  small  but  continuous  thrust.  Due  to  the 
geometry  of  an  EDT — we  will  assume  this  to 
consist  of  two  small  subsatellites  at  either  end 

Figure  1  -  An  EDT  uses  a  conductive  wire  to  interact 
J.  ^  ,  ,  ,  ,  .  with  the  geomagnetic  field.  Credit  Naval  Research  Labs 

of  the  long  conductive  tether — the  system  is 

dynamically  stable  in  a  nadir  pointing  orientation  due  to  gravity  gradient  effects.  The  nadir 
pointing  orientation  allows  thrust  generation  but  the  system  is  constrained  to  thrusting  in  the 
direction  of  the  cross  product  of  the  B  field  and  the  tether  orientation,  L.  This  thrust  direction 
constraint  combined  with  the  extremely  low  thrust  magnitude  pose  a  significant  challenge  in 


controlling  the  spacecraft.  Despite  the  system  complexity,  successful  control  of  an  EDT  has  been 
accomplished  in  simulations,  with  a  wide  range  of  orbital  maneuvers  being  achievable. 

While  an  inert  tether  is  dynamically  stable  in  its  nadir  pointing  orientation,  it  will  oscillate,  or 
librate,  both  in  and  out  of  the  orbital  plane  around  its  equilibrium  point.  The  motion  is  similar  to 
a  swinging  pendulum.  These  librations  are  due  to  the  competing  torques  from  gravity  gradient, 
atmospheric  drag,  and  if  the  tether  is  electrodynamic,  the  Lorentz  force.  Furthermore  Palaez  et  al. 
[7]  have  demonstrated  that  a  tether  with  a  constant  direct  current  will  eventually  become 
unstable.  Controlling  the  Lorentz  force  by  controlling  the  current  has  been  shown  to  limit  the 
amplitude  of  the  libration  to  within  acceptable  limits,  as  shown  by  Lanoix  et  al.  [8]  and  Hoyt  [9]. 
Williams  [10]  demonstrated  combined  libration  control  and  maneuvering  control.  He  derived  a 
numerical  method  for  determining  the  optimal  control  of  an  EDT  to  maneuver  between  two 
desired  orbits  and  also  demonstrated  that  the  natural  librations  can  be  exploited  to  achieve 
greater  thrust  by  taking  advantage  of  the  changing  tether  orientation.  This  method  works  well  for 
short  duration  maneuvers  as  he  uses  the  instantaneous  non-linear  equations  of  motion,  however; 
it  requires  significant  computation  time  as  a  result.  For  longer  scale  maneuvers,  on  the  order  of 
weeks  or  months  instead  of  hours,  this  method  is  both  infeasible  computationally  and  subject  to 
significant  accumulated  round  off  error.  As  a  consequence  of  the  low  thrust,  the  instantaneous 
state  will  only  vary  slightly  from  the  average  state,  which  slowly  changes  over  many  orbit 
revolutions.  By  applying  the  method  of  averaging  to  the  equations  of  motion,  Tragesser  and  San 
[11]  were  able  to  generate  sub-optimal  solutions  to  long  duration  maneuvers  by  using  a  periodic 
controller  to  change  the  average  state  of  the  system.  Their  work  assumes  a  non-librating,  nadir¬ 
pointing  tether  and  no  attempt  was  made  to  optimize  the  maneuvers.  Stevens  [12]  took  this 
method  of  averaging  and  combined  it  with  optimal  control  theory  to  generate  optimal  control 
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solutions  to  a  long  duration  tether  orbital  maneuver  while  considering  libration.  This  approach 
utilizes  a  mathematical  model  of  the  system  (i.e.  dynamic  equations  of  motion)  and  with 
optimization  software  numerically  determines  the  minimum  cost  (i.e.  time,  or  fuel  consumption) 
trajectory  to  achieve  the  desired  end  state.  The  validity  of  this  method  is  subject  to  the  quality  of 
the  mathematical  models  it  uses.  A  solution  generated  with  no  consideration  of  atmospheric  drag 
in  the  model  will  not  produce  reliable  performance  in  practice.  Likewise,  a  solution  generated 
with  a  poor  model  of  atmospheric  density  will  not  produce  reliable  results  when  implemented  in 
orbit.  Our  ability  to  accurately  predict  the  atmospheric  density  in  LEO  is  very  poor.  There  are 
highly  sophisticated  statistical  models  that  are  very  useful  in  satellite  tracking  and  orbit 
determination  but  these  models  do  not  reliably  predict  the  atmospheric  density  several  weeks 
ahead  of  time.  As  such,  the  optimal  control  solution  generated  using  these  models  will  not 
reliably  maneuver  the  spacecraft  to  the  desired  orbit.  The  purpose  of  this  paper  is  to  extend  the 
work  of  Stevens  by  implementing  feedback  control  to  the  optimal  control  maneuver  in  order  to 
achieve  reliable  performance  despite  modeling  inaccuracies. 

This  goal  will  be  achieved  in  four  steps  as  follows:  (also  visually  represented  in  Figure  3): 

1)  Derive  the  equations  of  motion  for  a  tether  system  from  the  Gaussian  variation  of 
parameters  equations. 

2)  Numerically  solve  the  non-linear  optimal  control  problem  (NLOCP)  for  the  whole 
maneuver. 

3)  Linearize  the  averaged  equations  of  motion  and  determine  the  linear  quadratic  regulator 
control  law  to  accommodate  a  dynamic  model  with  a  time  varying  atmosphere. 

4)  Test  the  feedback  controller  within  a  numerical  simulator  to  demonstrate  reliable 


performance. 
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To  produce  a  dynamic  model  we  will  assume: 

•  Near  circular  orbits  (e«l). 

•  Negligible  sensor  error  in  feedback  (complete  state  feedback). 

•  Perfect  knowledge  of  true  anomaly  (true  anomaly  as  the  independent  time  variable). 

•  Constant  nadir  pointing  tether  orientation 

•  Periodic  current  control  to  match  the  periodic  orbit  environment 

•  Non-tilted  dipole  model  of  geomagnetic  field 


Figure  3  -  Block  diagram  depiction  of  inner/outer  loop  structure 

II.  Derivation  of  Equations  of  Motion 

The  main  force  driving  the  equations  of  motion  is  central  body  gravitation.  We  will  also 
consider  the  influence  of  the  next  largest  force  in  the  low  earth  orbit  environment  which  is 
atmospheric  drag.  To  describe  an  orbit  we  require  six  pieces  of  information  that  comprise  the 
state  vector  or  element  set.  The  orbit  can  be  defined  as  a  state  vector  in  the  form  of  position  and 
velocity  vectors  or  as  an  element  set  in  the  form  of  a  scalar  magnitude  and  angular  positions. 
Position  and  velocity  vector  representations  are  not  the  most  efficient  method  for  describing 
orbits  as  the  state  values  will  constantly  be  changing  due  to  periodic  motion.  Instead  we  choose 
an  element  set  representation,  in  particular,  the  classical  orbital  elements  (we  in  fact  use  a  partial- 
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equinoctial  element  set,  a  distinction  to  be  described  later).  In  an  unperturbed  orbit  that  is  only 
subject  to  gravitation  influence,  the  classical  orbital  elements  (COEs)  do  not  vary  with  the 
exception  of  true  anomaly.  A  visual  representation  of  the  six  orbital  elements  is  seen  in  Figure  4 

[13]. 

PI 


Figure  4  -  The  six  classical  orbital  elements,  a,  e,  i,  £1 ,  w,  v.  credit:  Brandir 
These  six  states  completely  and  uniquely  define  an  orbit.  The  semi-major  axis  (a)  defines 

the  size  of  the  orbit,  eccentricity  (e)  defines  the  elongation  of  the  orbit,  inclination  (i)  defines  the 

orientation  of  the  orbit  with  respect  to  the  equator,  argument  of  perigee  (co)  defines  where  the 

lowest  point  of  the  orbit  is  with  respect  to  the  line  of  nodes,  right  ascension  of  the  ascending 

node  (f2)  defines  the  location  of  the  ascending  equatorial  crossings  (line  of  nodes)  and  true 

anomaly  (v)  defines  where  the  satellite  is  in  the  orbit  with  respect  to  perigee.  They  have  been 

carefully  chosen  to  be  static  in  value  during  free  motion.  In  forced  or  perturbed  motion,  Vallado 

[14]  defines  these  variation-of-parameters  (VOP)  equations  which  describe  their  time  rates  of 

change  due  to  any  perturbing  accelerations  (the  Gaussian  form).  Beginning  with  equation  (1.1),  a 

full  derivation  of  the  dynamics  can  be  found  in  Appendix  A 
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e~  V  1-e2  2 


h 

+  ecosv 


COSV/* - - - /s 
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(1.1) 


dt  r 

To  avoid  confusion  with  the  semi-major  axis  (a)  or  with  forces  ( F)  we  notate  the  perturbing 


accelerations  as  a  lowercase/ 


Figure  5  -  The  RSW  frame.  Radial  (R),  In-track  (S),  Orbit  normal  (W). 

For  the  equation  (1.1)  the  accelerations  must  be  expressed  in  the  Local  Vertical  Local 

Horizontal  (LVLH)  frame,  sometimes  also  written  as  RSW.  These  axes  are  directed  in  the  radial 

(R),  in-track,  (S),  and  orbit  normal  (W)  directions. 

1.  Equinoctial  Elements  -  h  and  k 

Our  first  step  is  to  modify  our  element  set.  We  do  this  because  of  singularities  for  zero 
eccentricity  and  the  poor  behavior  of  co  and  e  for  very  small  eccentricities.  These  considerations 
are  very  relevant  because  a  working  assumption  is  that  the  EDT  orbit  is  at  any  time 
approximately  circular  (e~0).  We  replace  the  elements  at  and  e  with  two  elements  borrowed  from 
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the  full  equinoctial  set  described  by  Vallado  [15].  We 
choose  h  and  k,  defined  as,  h  =  e sinvand  k  =  ecosv  . 
The  classical  pair,  eccentricity  and  argument  of 
perigee,  define  a  vector  (magnitude  and  direction 
respectively)  that  points  towards  perigee  of  the  orbit, 
see  Figure  6.  For  very  small  e ,  this  vector  can  rapidly 
cross  over  the  origin  causing  both  e  and  co  to  rapidly 


Evolution  of  eccentricity  vector  over  time 


k(-) 


change  in  an  almost  discontinuous  manner  (see  Figure  Figure  6  -  As  the  eccentricity  vector  crosses 

near  or  over  the  origin,  the  behavior 

7  and  8).  For  closer  passes  with  the  origin,  the  becomes  increasingly  singular  while 

equinoctial  elements,  h  and  k  vary  linearly. 

behavior  would  cause  numerical  errors,  eventually  becoming  non-differentiable. 


Time  (s)  Time  (s) 

„.  „  ,  .  ..  .  ..  .  .  Figure  8  -  Discontinuous  co  around  the  origin 

Figure  7  -  Discontinuous  eccentricity  around  the  origin 

The  equinoctial  elements  by  contrast  represent  the  Cartesian  components  of  this  vector  which 
avoid  singularities  when  the  orbit  is  circular.  When  the  eccentricity  vector  passes  by  the  origin, 
the  equinoctial  elements  exhibit  no  discontinuities.  By  avoiding  these  singularities  and  rapidly 
changing  values,  the  well  behaved  equinoctial  elements  h  and  k  are  much  more  suitable  for 


accurate  numerical  analysis. 


dco/dt  (rad/s) 
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With  this  element  substitution  in  mind  we  can  now  express  the  rates  of  change  for  the  new 
elements.  As  a  further  reduction  we  apply  our  assumption  of  nearly  circular  orbits,  whereby  we 
ignore  eccentricities  0{e2)  and  smaller.  In  deriving  the  equations  of  motion  we  will  eventually 


wish  to  find  the  averaged  state  dynamics  which  will  require  integration  of  the  dynamic 
equations.  To  simplify  the  process  of  integration  we  eliminate  time  and  allow  true  anomaly,  the 
angular  position,  to  be  the  independent  variable.  Changing  the  independent  variable  supposes  the 
system  has  perfect  knowledge  of  the  true  anomaly,  which  we  take  as  a  simplifying  assumption. 
The  resultant  dynamics,  equation  (1.2),  describe  the  changes  in  the  element  set  due  to  any 


perturbing  accelerations  in  the  RSW  frame. 


da 

dv 

di 

dv 

(Kl 

dv 

dh 
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dk 
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Me  sm  i 
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—  cos(v  +  cd)  —  h  sin  2v  —  2k  cos2  V 

T 

’/*" 

2sin(v'  +  <y)  +  /isin2  v-5k  sin  ccosc  — 4/z  cos2  v 

fs 

-k  sin(v  +  <y)cot/ 

.fw. 

sin(v  +  co)-k  sin  2v  -  2h  cos2  v 

T 

’/*" 

2  cos(v  +  cd)  +  k  sin2  v  +  5  h  sin  v  cos  v  —  4  k  cos2  v 

fs 

/z  sinfv  +  <y)  cot ; 

.fw. 

(1.2) 


To  continue  we  must  now  define  the  perturbing  accelerations. 


2.  Perturbing  Accelerations 

There  are  two  major  sources  of  perturbing  accelerations  we  will  consider: 
•  Electrodynamic  forces 


Atmospheric  drag  forces 
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We  first  will  consider  the  electrodynamic  forces  on  the  tether.  Using  a  non-tilted  dipole  model 
of  the  earth,  equation  (1.3),  and  our  assumption  that  the  tether  remains  nadir  pointing  we  derive 
the  electromagnetic  perturbing  accelerations  on  the  tether  (1.4)  by  the  Lorentz  force  definition, 
F  =  /L  rB . 


B  = 


-2sin(v  +  0))  sin  i 
cos(v/H-w)sinl 
cos  i 


(1.3) 


With  the  non-tilted  dipole  we  express  the  perturbing  accelerations  as: 


=  (l  +  3ecosv) 

0 

fs 

-cos  i 

Jw  J 

tether 

ma 

cos(v  +  co)  sin  i 

The  acceleration  is  a  function  of  the  state  vector,  x  =  [a  i  Q.  h  k\  ,  true  anomaly  (v), 
and  the  current  through  the  tether,  I  (the  control).  The  remaining  parameters  are  constants:  tether 
length  (L)  ,  tether  system  mass  (m)  ,  magnetic  constant  ( ju0  =  Ati-\ 0“7  N  ■  A  2 )  ,  and  the  dipole 

moment  of  the  earth  ( ym  =  4-1 022  A  •  m2 ) . 


Now  we  consider  atmospheric  drag  forces.  Using  the  standard  definition  of  drag 

I  C  A 

acceleration — F,  lm  =  — B  pV2\  ,  where  B*  =  — - — is  the  inverse  ballistic  coefficient — and 
g/  2  m 

the  flight  path  angle  to  define  the  velocity  vector  in  terms  of  orbital  elements,  we  can  describe 
the  atmospheric  drag  as  a  function  of  the  orbital  elements.  Again,  ignoring  0{e2)  and  smaller 
terms  we  find  the  accelerations  due  to  drag  are: 


fR 

1  n*  M 
=  --B  p 

2  a 

esinv 

fs 

l  +  2ecosv 

Jw. 

Drag 

0 

(1.5) 
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With  these  perturbing  accelerations  defined  we  write  out  fully  the  system  dynamics  as  rates 
with  respect  to  the  fast  time  variable  v  and  in  terms  of  the  control  current  I.  We  capture  the 

constant  values  in  Cm  =  . 

mnm 


—  =  2 Cm  cos  z‘(l+2e cosv)-7 
dv 

di  -C  .  .  2 

—  =  — -sin  z  cos  (v  +  co)  I 
dv  a 

dkl  -Cm  .  .  .  . 

— ^sin(v  +  a> )  cos(v  +  0))- 1 


dv  a 

dh  _  C„ 
dv  a 

dk  _  C„ 
dv  a 


-cost 


3 h  2 h  2k  . 

- 1 - cos  v  H - sint/  + 

2  e  e 


3k  2k  2 h  . 

—  - cosv - sint/  + 

2  e  e 


fh  hk 2A 

- 1 - V 


2  e 


cos  2v  + 


k  kh 


J 
2  y 


k  k(k2-h2Y 

2 +  2e2 

v  J 


sin  2v 


■I 


V  2 


cos  2v - 


h  h(k2-h2)^ 

2+  2e2 


sin  2v) 


(1.6) 


We  can  also  write  the  dynamics  in  a  much  more  compact  form,  which  we  will  use  for  the 
remainder  of  this  paper  where  the  state  is  denoted  as  x  and  Ai  and  Bi  are  matrices  defined  in 
Appendix  A. 


dx 

—  =  A1(x,v)  +  B1(x,v)7 
dv 


(1.7) 


II.  Method  of  Averaging 


So  far  the  state  dynamics  are  expressed  as  rates  with  respect  to  a  fast  time  variable,  the  true 
anomaly.  We  call  this  “fast”  time  because  it  describes  variations  that  occur  within  a  single  orbit 
on  a  short  time  scale.  Numerical  modeling  with  a  fast  independent  variable  is  problematic  for 
very  long  duration  maneuvers.  To  accurately  model  each  orbit  and  have  a  high  enough 
“resolution”  to  the  orbit  many  data  points  need  to  be  solved  for  within  each  orbit.  Maneuvers 
consisting  of  thousands  of  orbits  then  require  tens  to  hundreds  of  thousands  of  node  points.  Over 
such  long  time  spans,  numerical  round-off  errors  compound  and  can  become  significant. 
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Additionally,  the  computation  time  for  solving  the  optimal  control  problem  with  thousands  of 
nodes  is  infeasibly  large.  Because  the  orbit  perturbations  are  small,  orbit  variations  are  small 
over  short  time  spans  and  only  exhibit  observable  secular  changes  over  long  time  spans.  It  would 
be  computationally  advantageous  to  not  model  the  periodic  variations  at  all,  and  focus  solely  on 
the  secular  changes  in  the  state  as  in  Figure  7. 


Figure  7  -  Averaging  the  dynamics  over  N  integral  orbits. 


This  can  be  achieved  through  the  method  of  averaging.  The  secular  state  changes  occur 
slowly  enough  that  we  can  consider  them  constant  over  some  small,  integral  number  of  orbits,  N. 
Averaging  the  dynamics  over  that  time  span  has  the  effect  of  removing  the  periodic  variations. 


dx  1  c  dx. 

- = -  — dv 

dT  2 kN  {  dv 


(1-8) 


The  independent  variable  in  the  resulting  equations  is  no  longer  a  fast  time  variable  (v),  but  is 
now  a  slow  time  variable  (T).  The  slow  Equation  (1.8)  now  describes  the  secular  changes  in  state 
values.  These  averaged  states  have  no  real  meaning  at  any  single  instant  in  time,  but  serve  to 


capture  the  changes  in  the  system  state  over  time. 
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To  compute  the  integral  in  (1.8)  requires 


knowledge  of  the  variable  current,  1.  Looking  at 


equation  (1.2)  we  see  the  accelerations,  /,  are 


multiplied  by  periodic  functions  that  oscillate  at  the 


orbital  frequency  and  twice  the  orb  ital  frequency  Figure  8  -  The  square  wave  analog  to  a 

Fourier  series. 


(cosv  or  sinv  and  cos2v  or  sin2v).  To  get  the 

maximum  change  (i.e.  never  work  against  ourselves)  we  need  to  “push”  in  phase  with  the  natural 
dynamics  just  like  pushing  a  swing  in  sync  with  its  motion.  To  do  this  we  describe  the  current  as 
a  square  wave  analog  to  a  Fourier  series.  The  square  wave  aspect  arises  by  using  the  sign  of 
cosine  and  sine  functions  as  opposed  to  the  cosine  and  sine  function  themselves. 


where  u  and  *P  are, 


1  =  1 . ul'v  =  ry'u 


(1.9) 


uf=[mj  u2  m3  ua  u5\  (1.10) 

*Fr=sign([l  cos(v)  sin(v)  cos(2v)  sin(2v)])  (1.11) 

and  Im  is  the  maximum  RMS  current.  The  control  coefficients,  ui  through  U5,  are  bounded  as: 

-1  <  u  j  <  1 

ie[2,5],  -\[l<uj  <V2 

The  difference  between  Ui  and  U2  through  U5  is  because  the  ui  term  applies  to  the  DC  term 
while  the  m  through  U5  terms  apply  to  the  AC  sinusoidal  terms.  This  satisfies  the  RMS  current 
constraint, 


I 


2 

RMS 


=  U ,  H - 


V..2 


</ 


2 

RMS ,  Limit 


(1.12) 


19 


Again  note  the  sign  function  in  equation(l.ll),  this  is  a  departure  from  previous  work  in  the 
field  by  Stevens  [16].  In  his  work  Stevens  assumed  a  sinusoidal  periodic  controller.  We  have 

found  a  -27%  improvement,  a  factor  of  ["— j ,  by  using  square  wave  periodic  controller.  Visually 

the  improvement  is  seen  as  the  increased  area  under  the  square  wave  and  above  the  sine  or 
cosine  wave  (Figure  8).  As  the  cosine  or  sine  function  approaches  zero  the  current  would 
decrease  in  a  sinusoidal  controller  while  this  square  wave  controller  keeps  the  control  input 
(current)  high  until  the  natural  switching  point  of  the  dynamics  at  which  point  the  natural 
tendency  of  the  dynamics  reverses  and  so  the  influence  of  the  coefficient  switches  from  positive 
to  negative  (as  sinv  passes  from  positive  to  negative  the  effect  of  its  u  coefficient  switches  sign). 

With  the  system  dynamics  fully  described  by  (1.2),  (1.4),  and  (1.5),  we  now  undertake  the 
method  of  averaging.  The  matrices  found  below  are  explicitly  defined  in  Appendix  B.  By 
combining  those  three  equations  and  rearranging  we  can  express  the  state  dynamics  in  a  matrix 
form. 

—  =  A1(x,v)  +  B](x,v)'Pr(v)u  (1.13) 

dv 

Vector  Ai  represents  the  changes  due  to  uncontrolled  perturbations,  i.e.  drag.  Vector  Bi 
represents  the  changes  due  to  electrodynamic  forces  caused  by  some  current  I  in  the  tether.  It  is 
possible  to  extract  the  fast  time  dependence  from  Ai  and  Bi,  resulting  in  A'(x)<D(v)  and 
B  '(x)O(v)  .  This  simplifies  the  integration  by  only  requiring  integration  of  periodic  functions  of 
the  fast  time  variable.  Most  of  the  terms  in  the  product  average  out  to  zero.  Only  non-periodic 
secular  terms  survive  the  integration. 

The  resulting  averaged  dynamic  equations  of  motion  are  as  follows: 


—  =A2(x)  +  B,(x)u 
dT 


0.14) 
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where 

2  TIN  2  kN 

A,  = — - — A’,  [  <S>{v)dv  and  B,  = — i — B’,  f  ®(v)Tr(v)<iv-u 
2  In  N  1  J  2  2nN  1  J 

o  o 

Complete  expressions  for  the  preceding  equations  can  be  found  in  Appendix  B.  Equation 
(1.14)  will  serve  as  the  basis  for  the  remaining  two  steps. 

III.  Numerically  Solve  Nonlinear  Optimal  Control  Problem 

The  optimal  control  problem  (OCP)  is  to  find  the  control  vector  that  will  maneuver  from  a 
starting  point  to  a  final  end  point  while  minimizing  some  cost  of  the  problem.  Common  OCP’s 
are  minimum  fuel  maneuvers;  however  as  EDTs  are  not  fuel  limited  we  seek  minimum  time 
orbit  transfers.  To  formally  pose  the  OCP  we  must  find  the  control  function  that  will: 


Minimize  the  cost  function, 

II 

^5 

(2.1) 

subject  to  the  dynamic  constraints 

X  =  f  (x.u) , 

(2.2) 

boundary  constraints, 

o 

II 

9- 

(2.3) 

and  path  constraints, 

g(x,u)  =  0. 

(2.4) 

We  then  construct  the  Hamiltonian  as 

H  =XT  f 

(2.5) 

and  the  auxiliary  function 

<t>  =  (j)  +  vrcp 

(2.6) 

where  X  is  the  costate  vector  and  v  is  a  Lagrange  multiplier.  The  necessary  conditions  for 


optimality  are  then  the  adjoint  equations 

l--dH 

ax 

(2.7) 

and  the  control  optimality  condition 

^  =  0 

du 

(2.8) 

subject  to  the  transversality  conditions  A(tf )  =  -  —  \t=t 


(2.9) 
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and  +  lf=(/  =  0  (2.10) 

From  these  necessary  conditions  we  have  a  two  point  boundary  value  problem.  Conceptually 
we  have  a  set  of  coupled  partial  differential  equations  consisting  of  state  dynamics  (equation 
(2.2))  and  costate  dynamics  (equation  (2.7))  with  various  initial  or  final  conditions.  There  are 
many  methods  for  solving  such  systems  numerically.  We  implement  a  program  called  DIDO  to 
find  a  candidate  optimal  control  function.  The  reason  we  chose  this  software  is  that  it  provides 
quick  and  accurate  solutions  without  a  first  guess  of  the  solution.  Additionally  the  software  is 
easy  to  use  and  it  does  not  require  the  user  to  fully  understand  pseudospectral  (PS)  methods  (the 
method  used  internally  to  find  a  solution).  The  user  simply  has  to  write  the  OCP  much  as  they 
would  on  paper  (or  as  above).  The  PS  method  does  not  use  propagation  and  thus  it  is  not 
susceptible  to  the  round-off  errors  associated  with  other  methods.  For  a  full  description  of  PS 
methods  and  their  benefits  see  Refs  [17],  [18]. 

IV.  Linear  Control  Loop 

The  optimal  control  solution  is  a  pure  open-loop  control  and  assumes  a  perfect  dynamic 
model  of  the  system.  In  reality  however,  poor  density  and  geomagnetic  field  models  can  lead  to 
significant  errors.  To  combat  these  inaccuracies  and  ensure  reliable  performance  we  introduce 
linear  feedback  as  shown  in  Figure  3.  Due  to  the  exceptionally  small  thrust  produced  by  an  EDT, 
maneuvers  happen  over  very  long  time  scales.  The  states  change  slowly  and  can  be  considered  to 
be  constant  within  a  single  orbit  period.  With  this  assumption  we  linearize  the  dynamics  from 
(1.14)  with  the  intent  of  generating  a  linear-quadratic  regulator  (LQR)  control  law.  The 
linearized  dynamics  we  generate  here  will  have  state  dependent  gains.  To  maintain  this  as  a 
linear  time  invariant  system  we  update  the  gains  once  per  orbit.  The  computation  time  for  this 
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method  is  short  compared  to  the  hold  time,  with  at  most  several  seconds  for  the  computation  and 
90  minutes  of  orbit  at  minimum.  This  avoids  treating  the  system  as  linear  time  varying  by  doing 
a  zero  order  hold — computing  the  gains  at  the  beginning  of  the  orbit  with  the  states  at  that  time 
and  not  updating  until  the  next  orbit. 

Using  the  following  procedure  (fully  developed  in  Appendix  C)  we  linearized  the  averaged 
equations  of  motion. 

x  =  f(x,u)  (3.1) 

Making  a  first  order  approximation  of  we  say  that 

f  (x  +  8x,  u  +  Su)  =  f  (x,  u)  +  f  tSx  +  f„Su  (3.2) 

We  can  reduce  this  to  the  following 

Sx  =  frSx  +  f„Su  (3.3) 

The  subscripts  on  f  denote  the  Jacobian  of  f  with  regards  to  the  state  or  control  vector,  as 
appropriate.  For  the  linear  system.  Equation  (3.3),  with  the  following  quadratic  cost  function 

'/ 

/  =  |  (xrQx  +  i/Rujc/f  (3.4) 

The  feedback  control  law  that  minimizes  the  cost  is 

8u  =  -K8x  (3.5) 

where  K  is  the  solution  the  linear  quadratic  regulator  problem  and  is  given  by 

K  =  R  VP  (3.6) 

and  P  is  the  solution  to  the  continuous  time  Riccati  differential  equation. 

As  we  are  assuming  no  sensor  feedback  we  are  able  to  define  a  state  space  system  with  A  =  f  r , 

B  =  f„ .  Upon  proper  definition  of  weighting  matrices,  Q  and  R  this  is  now  sufficient  information 
to  solve  for  the  LQR  gain  matrix,  K.  We  generate  a  solution  using  the  built-in  LQR(...)  function 


in  MATLAB  7.10.0  (R2010a). 
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The  candidate  solution  output  from  DIDO  consists  of  averaged  states  thus  a  simple 
interpolation  of  those  would  not  give  us  an  appropriate  reference  trajectory.  The  instantaneous 
trajectory  is  known  to  have  periodic  oscillations  (the  entire  reason  for  averaging  to  begin  with), 
and  so  to  generate  the  reference  trajectory  we  propagate  forwards  one  orbit  using  the  state 
predicted  by  DIDO  as  the  initial  state.  This  method  of  on-the-fly  reference  generation  avoids  the 
numerical  issues  that  arise  if  the  entire  reference  is  propagated  from  the  beginning.  Numerical 
propagation  over  a  significant  number  of  time  steps  accumulates  round-off  errors.  For  maneuvers 
lasting  100  orbits  and  each  orbit  consisting  of  roughly  5-10,000  time  steps  there  is  considerable 
room  for  accumulated  errors.  The  problem  is  only  compounded  as  the  maneuver  duration  is 
increased.  In  our  method  the  propagating  time  span  is  at  most  one  orbit  period  thus  there  is  no 
significant  accumulated  round  off  error. 

V.  Results  &  Analysis 

To  test  the  control  law  we  must  simulate  the  system  numerically  because  we  cannot  recreate 
the  EDT  dynamics  physically  on  the  ground.  Using  MATLAB  we  created  a  simulator,  the 
structure  of  this  process  is  captured  in  Figure  9.  This  is  essentially  identical  to  the  block  diagram 
presented  earlier  (Figure  3),  the  main  difference  being  the  “Plant”  is  now  simply  another 
numerical  propagator  with  modeling  error  introduced  rather  than  a  physical  system. 
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Figure  9  -  Flowchart  of  Control  Law  testing.  Generate  the  control  law  and 
also  generate  the  corrected  reference.  Iterate  through  the  Orbit  Loop  until 
one  full  orbit  has  passed  then  restart  the  process  for  the  next  orbit. 


During  numerical  testing  (using  the  above  scheme)  we  sought  to  characterize  the  performance 
of  a  linear  controller  when  atmospheric  density  model  errors  were  introduced.  To  implement 
these  errors  we  applied  a  lognormal  distribution  to  the  density.  Measured  data  for  altitudes 
ranging  from  120km  up  to  750km  have  proven  to  be  consistent  with  a  lognormal  distribution 
rather  than  the  more  common  normal  distribution,  [19],  [20].  The  lognormal  distribution  is 
appropriately  named  because  the  log  of  the  values  is  normally  distributed.  Lognormally 
distributed  values,  X,  can  be  generated  from  a  set  of  normally  distributed  numbers,  Z,  by 

X  =  eZa+M  (4.1) 

where  a  is  the  standard  deviation  of  the  logarithmic  values,  ln(Y) ,  and//  is  their  mean.  The  mean 

is  obtained  in  a  table  look-up  of  the  1976  Standard  Atmosphere  and  the  standard  deviation  is 
fixed  at  1%  of  the  mean.  For  this  example  a  mean  of  10"  kg  m"  was  chosen  purely  for 


illustrative  purposes. 
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Lognormal  probability  density  function  of  p 


Figure  10  -  The  probability  density  function  for  density  is  a  log-normal  distribution. 


Considering  the  linear  system,  the  “noise”  due  to  density  variations  is  not  white  however;  the 
values  do  have  a  mean  of  zero  because  the  lognormal  values  only  ever  occur  with  a  -1  attached 
to  it,  shifting  them  over  0.  Comparing  this  shifted  lognormal  distribution  with  a  Gaussian 
distribution  for  low  values  of  standard  deviation  (<0.1)  we  find  that  the  two  distributions  are 
markedly  similar,  with  only  a  slight  skew  to  the  right  in  the  lognormal  distribution.  Under  this 
premise,  linear  quadratic  control  is  still  valid,  though  not  “optimal”  because  of  the  non-white 
noise.  Using  this  method  we  set  about  to  demonstrate  the  performance  under  the  following 
conditions. 

Maximum  altitude  raise  in  a  fixed  time  span: 

Initial  Altitude  :  250km 
Initial  Inclination  :  1  degree 
Initial  Eccentricity  :  .05 1 
Duration  :  20  orbits 
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o 

'ro 

E 

E 

CD 

CO 


8  10  12 
Time  (orbits) 


Figure  11  -  A  comparison  of  controlled  and  uncontrolled  responses  during  an  orbit  raise  maneuver. 


The  “controlled”  scenario  is  the  scenario  where  we  apply  linear  feedback  control  to  maintain 
the  reference  state.  The  “uncontrolled”  scenario  is  the  scenario  where  we  apply  no  feedback  and 
only  have  the  reference  optimal  control  as  input.  The  controlled  system  remains  bounded  around 
zero  error  while  the  uncontrolled  system  has  a  state  error  that  grows  unbounded.  After  20  orbits, 
which  is  roughly  one  and  a  half  days,  the  state  error  is  already  approaching  1.5  km. 

It  is  interesting  to  note  the  bias  in  the  error  towards  positive  altitude  errors.  We  can  explain 
this  by  again  looking  at  the  density  distribution.  The  lognormal  distribution  has  a  slight  skew  to 
the  right,  meaning  slightly  higher  densities  more  often.  Higher  densities  would  create  a  larger 
drag  force,  which  would  tend  to  reduce  the  orbit  altitude.  The  error  is  computed  as  the  difference 
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between  the  reference  and  the  measured  state,  so  for  altitudes  that  are  decreasing  faster  than 
expected,  the  state  error  tends  towards  the  positive.  This  causes  a  steadily  increasing  state  error 
in  the  uncontrolled  response  (the  error  is  “noisy”  but  it  is  overlaid  onto  a  roughly  linearly 
increasing  component). In  the  controlled  scenario  we  see  a  steady  state  bias  towards  the  positive. 
The  bias  is  again  due  to  the  non-white  (lognormal)  distribution  of  the  atmospheric  density  noise. 
The  LQR  method  only  generates  the  optimal  gain  matrix  for  white  noise  systems  so  the 
emergence  of  a  steady  state  error  is  due  applying  LQR  to  lognormal  noise. 

VI.  Conclusion 

We  have  successfully  applied  linear  feedback  control  principles  to  a  non-linear  optimal 
control  problem,  specifically  electrodynamic  tether  maneuvers.  We  have  demonstrated  a  flaw 
with  existing  non-linear  optimal  control  methodologies  due  to  model  uncertainties.  Under  solely 
the  influence  of  the  “optimal  control”  when  the  atmospheric  density  varies  lognormally,  there  is 
introduced  a  significant  state  error.  After  applying  linear  feedback  control  via  a  linear  quadratic 
regulator,  the  state  errors  remain  bounded.  This  method  allows  for  reliable  and  predictable,  if 
perhaps  slightly  sub-optimal,  performance  of  an  electrodynamic  tether. 

There  are  also  many  ways  to  extend  the  work  completed  here.  A  more  sophisticated 
geomagnetic  field  model  or  one  with  model  uncertainties  can  be  implemented.  Other  smaller 
perturbations  can  be  considered  such  as  solar  radiation  pressure  and  non-spherical  earth 
dynamics.  Additionally,  linear  feedback  control  methods  other  than  linear  quadratic  regulators 
can  be  investigated  for  better  performance  and  reduced  stead  state  error  when  subject  to 
lognormal  noise.  Such  improvements  will  further  enhance  the  reliability  this  method  provides  to 
precision  maneuvering  of  electrodynamic  tethers. 
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Appendix  A 


Derivation  of  Equations  of  Motion 

All  equations  references  are  from  “Fundamentals  of  Astrodynamics  and  Applications”  by  David 
A.  Vallado.  We  begin  with  the  Gaussian  VOP  equations  from  Vallado. 
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We  now  undergo  a  process  of  applying  several  transformations.  These  are: 

•  Replacement  of  e,  co  with  h,  k  (partial  equinoctial  set). 

o  h  =  e  sin(a>)  +  e  cos(a>) 
o  k  =  ecos(a>)  -  esin(&>) 

•  Nearly  circular  orbits  (e«l). 

o  This  allows  us  to  simplify  many  expressions  containing  0(e2)  terms.  Also 
terms  dividing  by  (1  -  e) become  multiplied  by  (1  +  e) etc. . . 

•  Replacement  of  time  as  independent  variable. 

o  Multiply  each  equation  of  the  system  dynamics  by  —  =  —  (l-2ecosv) 

dv  jU 

This  leads  us  to  the  reduced  dynamics  in  the  partial  equinoctial  set. 


32 


■y-  =  —  [esm(v)fR  +  (l-ecos(V))/x] 
dv  Me, 


di  a~ 
dv  Me 


(1  -  3e  cos(v))  cos(y  +  co)fw 


dv  Me  sin(i) 


2 

-cos (y  +  cd)  —  h  sin(2v)  -  2k  cos 2 (v) 

T 

"/*" 

Cl 

2  sin(v  +  co)  +  h  sin2  (v)  -  5k  sinfv)  cos(v)  -  4/z  cos2  (v) 

fs 

M® 

-k  sin(v  +  <y)cot(z) 

Jw_ 

2 

sin(v  +  co)  —  k  sin(2v)  -  2/z  cos2  (v) 

■T 

7/ 

Cl 

2cos(v  +  co)  +  k  sin2(v)  +  5/zsin(v,)cos(v/)  -4Tcos2(v) 

fs 

h  sinfv  +  co)  cot(z) 

Jw. 

The  accelerations  due  to  electrodynamic  interactions  are  a  result  of  the  Lorentz  force, 
F  =  TLaB  where  L  is  of  magnitude  L  in  the  Radial  direction  and  B  is  the  local  geomagnetic  field 


vector.  We  used  a  non-tilted  dipole  model,  B 


r,„M0 

3 


-2  sin  (v  +  co)  sin(/) 
cos(v  +  <y)  sin(z) 
cos(z) 


(N-A  'm1),  though  any 


model  could  be  used.  We  need  accelerations  rather  than  forces,  so  we  write  /  =F/m  to  express 
the  accelerations  in  the  LVLH  frame. 


f tether  =  (l+  3g  COS(V)) 


0 

—cos  (z) 

cos(v  +  co)  sin( ;) 


(m-s"2) 


To  characterize  the  accelerations  due  to  drag  we  first  express  the  force  as 


Fdrag  =  ~^CDApVn\  .  We  can  determine  the  velocity  at  any  point  via 


V  = 


tL 

a 


The  direction  of  the  velocity  vector  is  defined  by  the  flight  path  angle,  </>fpa .  This  is  the  angle 


measured  from  the  local  horizontal  to  the  velocity  vector. 
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V  = 
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(4.2) 


It  is  thus  unnecessary  to  define  (j)fpa  exactly  as  it  is  much  easier  to  define  the  sine  and  cosine  of 
this  angle.  From  Vallado  we  can  say: 
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From  this  we  can  fully  express  the  drag  acceleration.  Before  we  do  this,  we  introduce  a  new 

variable,  B  .  We  define  this  “inverse  ballistic  coefficient”  as  B  =  — ^  .  Finally  we  express  the 

m 


drag  accelerations  as 
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Putting  together  the  tether  and  drag  forces  with  the  general  system  dynamics  we  arrive  at  the 
following  equations, 
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Where  C  =  .  This  is  a  very  cumbersome  set  of  equations  to  work  with.  We  can  notice 

many  common  terms,  notably,  sine  and  cosine  of  v  or  2v.  If  we  extract  these  terms  and  write 
these  five  equations  as  matrices,  we  arrive  at  the  following. 
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Which  we  can  reduce  compactly  to  the  following  matrix  form  which  is  the  basis  of  the 
remaining  analysis, 

dx 

—  =  A1(x,v)  +  B1(x,v)7 

dv 
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Appendix  B 


Method  of  Averaging  -  Matrices 


The  purpose  of  this  appendix  is  to  define  the  matrices  found  when  applying  the  method  of 
averaging.  Those  matrices  would  be  A  and  B  and  their  many  variations,  and  the  product. 
Tr=sgn([l  cos(v)  sin(v)  cos(2v)  sin(2v)]) 


<I>r  =  [l  cos(v)  sin(v)  cos(2v)  sin(2v)] 


a  ae  0  0  0 
0  0  0  0  0 


A,j(x) 


-B  pa 


0 

h 

2 

k 


0 

h 

e 

k 


12  e 

A1(x,v)  =  A1'(x)0(v) 


0  0  0 


k 

e 


-h  -k 


e 


2  cos(i) 

— —  sinO') 

2  a 

4e  cos  O') 

0 

0 

0 

0 

k2-h2  .  ... 

2  sin(i) 

Zae 

hk 

0 

0 

0 

ae2 

3  h 

2  h 

2k  1 

f  h  hk 2  ^ 

- cos(i) 

—cos  (0 

—  cos(i) 

—  “• - 2 

Jcos(i) 

2  a 

ae 

ae 

\2 a  ae  j 

3  k 

2k 

2  h 

•  k  kh 2 ' 

cos(i) 

—  cos  (l) 

- cos(i) 

- cosp) 

2  a 

ae 

ae 

yla  ae  ) 

B1(x,v)  =  B1,(x)0(v)Tr(v) 


0 

hk  . 

— -sin(0 
ae 

k2-h2 


2  ae2 


k  k(k2-h2) 


2  a 


2  ae 


h  +  h(k2-h2) 


2  a 


2  ae" 


COS  O') 


cos(i) 


Now  we  apply  the  method  of  averaging,  whereby  the  A’  and  B’  matrices  are  not  dependent  on 

T 

time,  and  thus  can  be  pulled  out  of  the  integral.  This  leaves  us  to  integrate  O,  and  OT  . 


And  those  integrals  are: 
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1 

2nN 


2k  N 

j  <S>dv 

o 


1 

0 

0 

0 

0 


2xN 


2kN 


LJt  IS 

J  <f>(v)x¥T  (v)dv 


Z7TA 

I 


1 


cos(v) 

sin(v) 


sgn(cos(y)) 

|cos(v)| 

sgn(cos(v/))sin(v/) 


2  ItN 


2tcN 


sgn(sin(v)) 
sgn  (5111(12))  cos(v) 

|sin(v)| 

cos(2v)  sgn(cos(v/))cos(2v/)  sgn  (sin (y))  cos (2v) 
sin(2v)  sgn(cos(v/))sin(2v/)  sgn(sin(v/))sin(2v/)  sgn  (cos(2v/))sin(2v/) 

1  0  0  0  0 

0  2jn  0  0  0 

0  0  2 jn  0  0 

0  0  0  2jn  0 

0  0  0  0  2jn 


sgn(cos(2v))  sgn(sin(2v)) 

sgn  (cos(2v))  cos(v)  sgn  (sin(2v))  cos(v) 

sgn  (cos(2v/))sin(v/)  sgn  (sin(2v/))sin(v/) 

|cos(2v)|  sgn  (sin(2v/))cos(2v/) 

|sin(2v)| 


\dv 


J  <t>(v)'PT (v)dv  = 


That  lets  us  define  A2  and  B2  as  the  following: 


2  JtN 


A,  =- 


2  kN 


2  JtN 


A \  j"  <f>(v)dv  and  B2  = - B’,  [  ®(v),Pr {v)dv-n 

J  2nN  J 


A,(x)  =  -B  pa 


a 

0 

0 

-h 

2 

1  , 


B2(x)  =  ■^2L 
7ta 


27Tacos(i )  8aecos(i) 


-^srsin(i) 


0 

k2-h 2 


sin(i) 


2  hk 

2 


3  Trh  ...  4  h  ...  4  k 

- cos(0  — cos(0  — cos(i) 

2  e  e 

3  Jtk  4  k  4  h 

- COS(i)  - cos(i) - COS(z) 

lee 


rh  2hk2^ 
r  k  2 kh1 


1 


2 hk  .  ... 
— —  sin(i) 


k2-h2 


cos  (i) 


cos(0  - 


A 

2 

e 

c2-/r)^ 

2 

e 

/ 

+  'A 

k2-h2) 

cos(/) 


cos(/) 
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Appendix  C 


Linearization 


Beginning  with  the  averaged  equations  of  motion 

dx 

—  =  A2(x)  +  B,(x)u 
dT 

we  apply  a  first  order  approximation  for  f  (x  +  5x,  u  +  5u )  as  the  following: 


f  (x  +  5x,  u  +  5u)  =  f  (x,  u)  +  f  t5x  +  fu  5u 


fx  and  fu  are  the  Jacobian  of  dx/dT  with  respect  to  the  state  vector,  x,  and  the  control  vector,  u. 

d  (A^  +B,u) 

The  second  term  is  simple  enough  to  compute  by  had  because  it  reduces  to  B3  = - - - =  B,  because 

du 

neither  A2  nor  B2  are  functions  of  u. 

The  first  term  however,  is  significantly  more  difficult.  Both  A2  and  B2  are  functions  of  the  state  vector 
and  thus  have  filled  Jacobians.  This  operation  was  done  with  the  help  of  a  symbolic  manipulator  to  ensure 
no  human  errors  were  introduced  in  the  transcription  of  the  equations  due  to  their  complexity.  Also,  even 
in  simplified  expressions,  these  matrices  are  too  large  to  display  properly  on  a  page  and  for  that  reason 


d  (A,  +B,u) 

will  not  shown  in  “equation  form”.  The  best  definition  of  A3  is  as  A,  = - =  fx .  The 

dx 


matrix  is 


presented  element-wise  in  the  form  used  in  the  computer  code. 
First  Row,  Columns  1-5 

-2*Bs*p*a, 

-2*Cm*sin  (i)  *ul-8*Cm/pi*  (hA2+kA2)  A  (1/2)  *sin(i)  *u2, 

0, 

8*Cm/pi/ (hA2+kA2) A (1/2) *cos (i) *u2*h, 

8*Cm/pi/ (hA2+kA2) A (1/2) *cos (i) *u2*k; 


Second  Row,  Columns  1-5 

ls*Cm/aA2*sin  (i)  *ul-Cm/pi/aA2*  (-kA2+hA2)  /  (hA2+kA2)  *sin(i)  *u4-2*Cm/pi/aA2*h*k/  (hA2+kA2)  *sin(i)  *u5, 
-l/2*Cm/a*cos (i) *ul+Cm/pi/a* (-kA2+hA2) / (hA2+kA2) *cos (i) *u4+2*Cm/pi/a*h*k/ (hA2+kA2) *cos (i) *u5. 


0, 
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2*Cm/pi/a*h/ (hA2+kA2 ) *sin ( i) *u4-2*Cm/pi/a* (- 

kA2+hA2) / (hA2+kA2) A2*sin ( i ) *u4*h+2* Cm/pi /a*k/ (hA2+kA2) *sin (i) *u5- 
4*Cm/pi/a*hA2*k/ (hA2+kA2) A2*sin ( i ) *u5, 

-2 * Cm/pi/ a* k/ (hA2+kA2 ) *sin(i) *u4-2*Cm/pi/a* (— 

kA2+hA2) / (hA2+kA2) A2*sin(i) *u4*k+2* Cm/pi /a*h/ (hA2+kA2) *sin(i) *u5- 
4*Cm/pi/a*h*kA2/ (hA2+kA2) A2*sin(i) *u5; 


Third  Row,  Columns  1-5 


2*Cm/pi/aA2*h*k/ (hA2+kA2) *u4-Cm/pi/aA2* (-kA2+hA2) / (hA2+kA2) *u5, 

0, 

0, 

-2*Cm/pi/a*k/ (hA2+kA2) *u4+4*Cm/pi/a*hA2*k/ (hA2+kA2) A2*u4+2*Cm/pi/a*h/ (hA2+kA2) *u5-2*Cm/pi/a* (— 
kA2+hA2 ) / (hA2+kA2 ) A2*u5*h, 

-2*Cm/pi/a*h/ (hA2+kA2) *u4+4*Cm/pi/a*h*kA2/ (hA2+kA2) A2*u4-2*Cm/pi/a*k/ (hA2+kA2) *u5-2*Cm/pi/a* (- 
kA2+hA2 ) / (hA2+kA2 ) A2*u5*k; 


Fourth  Row,  Columns  1-5 


-l/2*Bs*p*h-3/2*Cm/aA2*h*cos (i) *ul-4*Cm/pi/aA2*h/ (hA2+kA2) A (1/2) *cos (i) *u2- 
4*Cm/pi/aA2*k/ (hA2+kA2) A (1/2) *cos (i) *u3-Cm/pi/aA2* (h+2*h*kA2/ (hA2+kA2) ) *cos (i) *u4- 
Cm/pi/aA2* (k+k* (kA2-hA2) / (hA2+kA2) ) *cos (i) *u5, 

-3/2*Cm/a*h*sin (i) *ul-4*Cm/pi/a*h/ (hA2+kA2) A(l/2)*sin(i) *u2- 

4*Cm/pi/a*k/ (hA2+kA2) A (1/2) *sin (i) *u3-Cm/pi/a* (h+2*h*kA2/ (hA2+kA2) ) *sin (i) *u4-Cm/pi/a* (k+k* (kA2- 
hA2 ) / (hA2+kA2 ) ) *sin (i) *u5, 

0, 

-l/2*Bs*p*a+3/2*Cm/a*cos (i) *ul+4*Cm/pi/a/ (hA2+kA2) A (1/2) *cos (i) *u2- 
4*Cm/pi/a*hA2/ (hA2+kA2) A (3/2) *cos (i) *u2- 

4*Cm/pi/a*k/ (hA2+kA2) A (3/2) *cos (i) *u3*h+Cm/pi/a* (l+2*kA2/ (hA2+kA2) - 
4*hA2*kA2/ (hA2+kA2 ) A2) *cos (i) *u4+Cm/pi/a* (-2*k*h/ (hA2+kA2) -2*k* (kA2- 
hA2 ) / (hA2+kA2 ) A2*h) *cos (i) *u5, 

-4*Cm/pi/a*h/ (hA2+kA2) A (3/2) *cos (i) *u2*k+4*Cm/pi/a/ (hA2+kA2) A (1/2) *cos (i) *u3- 
4*Cm/pi/a*kA2/ (hA2+kA2) A (3/2) *cos (i) *u3+Cm/pi/a* (4*k*h/ (hA2+kA2) - 

4*h*kA3/ (hA2+kA2 ) A2) *cos (i) *u4+Cm/pi/a* (1+ (kA2-hA2) / (hA2+kA2) +2*kA2/ (hA2+kA2) -2*kA2* (kA2- 
hA2 ) /  (hA2+kA2 ) A2 ) *cos (i) *u5;  .  .  . 


Fifth  Row,  Columns  1-5 


-l/2*Bs*p*k-3/2*Cm/aA2*k*cos (i) *ul- 

4*Cm/pi/aA2*k/ (hA2+kA2) A (1/2) *cos (i) *u2+4*Cm/pi/aA2*h/ (hA2+kA2) A (1/2) *cos (i) *u3-Cm/pi/aA2* (k— 
2*k*hA2/ (hA2+kA2 ) ) *cos (i) *u4-Cm/pi/aA2* (-h-h* (kA2-hA2) / (hA2+kA2) ) *cos (i) *u5. 


-3/2*Cm/a*k*sin (i) *ul- 

4 * Cm/pi/ a* k/ (hA2+kA2) A (1/2) *sin (i) *u2+4*Cm/pi/a*h/ (hA2+kA2) A (1/2) *sin(i) *u3-Cm/pi/a* (k— 
2*k*hA2/ (hA2+kA2 ) ) *sin (i) *u4-Cm/pi/a* (-h-h* (kA2-hA2) / (hA2+kA2) ) *sin (i) *u5, 

0, 

-4*Cm/pi/a*h/ (hA2+kA2) A (3/2) *cos (i) *u2*k- 

4*Cm/pi/a/ (hA2+kA2) A (1/2) *cos (i) *u3+4*Cm/pi/a*hA2/ (hA2+kA2) A (3/2) *cos (i) *u3+Cm/pi/a* (- 

4*k*h/ (hA2+kA2) +4*k*hA3/ (hA2+kA2) A2) *cos (i) *u4+Cm/pi/a* (-1- (kA2- 

hA2 ) / (hA2+kA2) +2*hA2/ (hA2+kA2 ) +2*hA2* (kA2-hA2) / (hA2+kA2) A2) *cos (i) *u5, 

-l/2*Bs*p*a+3/2*Cm/a*cos (i) *ul+4*Cm/pi/a/ (hA2+kA2) A (1/2) *cos (i) *u2- 

4*Cm/pi/a*kA2/ (hA2+kA2) A (3/2) *cos (i) *u2+4*Cm/pi/a*k/ (hA2+kA2) A (3/2) *cos (i) *u3*h+Cm/pi/a* (1- 
2*hA2/ (hA2+kA2) +4*hA2*kA2/ (hA2+kA2) A2) *cos (i) *u4+Cm/pi/a* (-2*k*h/ (hA2+kA2) +2*k* (kA2- 
hA2 ) / (hA2+kA2 ) A2*h) *cos (i) *u5]  ; 


